Using dynamically scattered electrons for 3-dimensional potential reconstruction 
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Three-dimensional charge density maps computed by first-principles methods provide information 
about atom positions and the bonds between them, data which is particularly valuable when trying 
to understand the properties of point defects, dislocations, and interfaces. This letter presents a 
method by which 3-dimensional maps of the electrostatic potential, related to the charge density map 
by the Poisson equation, can be obtained experimentally at 1 A resolution or better. This method 
requires data acquired by holographic transmission electron microscopy (TEM) methods such as 
off-axis electron holography or focal series reconstruction for different directions of the incident 
electron beam. The reconstruction of the 3-dimensional electrostatic (and absorptive) potential is 
achieved by making use of changes in the dynamical scattering within the sample as the direction 
of the incident beam varies. 
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Extended defects such as dislocations, interfaces, in 
particular solid liquid interfaces, pose a particular chal- 
lenge to atomic scale computational methods. While 
their complexity normally requires super-cells too large 
for ab-initio methods, the long-range nature of forces de- 
termining their properties (e.g. strain fields, Coulomb 
and van der Waals [dispersion] forces) makes the con- 
struction of reliable yet efficient interatomic potentials 
required for molecular dynamics (MD) and related com- 
putational methods an extremely complicated task. Be- 
ing able to perform experiments which directly map the 
3-dimcnsional local electrostatic potential would provide 
a wealth of information without the need to do any sim- 
ulations at all and, if computations are needed to extract 
information not provided by the potential map or, in the 
case of MD, it would allow verification of interatomic 
potentials by direct comparison between experiment and 
simulation for the very (defect) structure under investiga- 
tion. For single crystals of small unit cell it has recently 
been demonstrated that it is possible to map the bond 
charge distribution by fitting Fourier coefficients of the 
crystal potential to convergent beam electron diffraction 
(CBED) data Q. 

TEM images, in the context of focal series reconstruc- 
tion (inline holography) over a large defocus range [lij ]. 
are very sensitive to relative atomic positions and varia- 
tions in the electrostatic potential. Although claiming to 
be able to provide data on the 3-dimcnsional distribution 
of atoms this is not true for destructive methods such as 
3D atom probe 0] because the atoms being " imaged" are 
not extracted from their original environment, but from 
the sample surface. 

The enormous advantages of TEM imaging come at 
the price of the lack of a general direct interpretability at 
atomic resolution. Although modern electron optics has 
been able to remove many of the artefacts caused by aber- 
rations of the imaging system, it cannot remove artefacts 
produced by the multiple (or dynamical) scattering pro- 



cess itself. This is the main reason why electron tomogra- 
phy has never been applied at atomic resolution, with the 
exception of extremely small nanostructures consisting of 
light atoms, for which the authors felt that they could ne- 
glect dynamical electron scattering events [3]. Attempts 
to correct for artefacts (i.e. features which cannot be in- 
terpreted directly in terms of atomic structure) in the exit 
face wave function produced by dynamical scattering ef- 
fects have, in the case of thin or non-crystalline samples, 
at most been able to reconstruct an average projected 
potential |, H, @, 0, S] • 

The method presented in this letter directly recon- 
structs the local complex scattering potential on a 3- 
dimensional grid, the real part of which is the electro- 
static potential, by separating the multiple scattering 
paths between potential voxels. 

The following discussion will, in order to keep the equa- 
tions readable, only consider the reconstruction of a 2- 
dimensional structure from a series of 1-dimensional im- 
ages. The extension to 3 dimensions is straight forward, 
requiring only small changes in the presented formalism. 

The proposed reconstruction method is based on the 
real-space variant 0] of the multislice algorithm [l(| for 
solving the relativistically corrected Schrodinger equa- 
tion describing the propagation of a fast electron through 
the specimen potential (see figure Q] for an illustration). 
In a first step the 2-dimensional scattering potential (re- 
member that this discussion can easily be extended to 3 
dimensions) is divided into a set of N discrete horizontal 
slices of thickness e (the optical axis of the microscope is 
assumed to be vertical, the fast beam electrons travelling 
down the microscope column) and the potential within a 
slice of index m is projected into a 1-dimensional layer of 
potential 
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The electron propagation is then approximated 
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FIG. 1: (color online) Diagram of real-space multislice algo- 
rithm. The scattering of the incident electron beam described 
by a plane wave is approximated by a finite number of scat- 
tering events at equidistant layers partitioning the sample in 
z-direction. As a consequence, the signal in the exit face wave 
function is non-local, including contributions from a number 
of scattering paths, the relative phase of which can be varied 



by changing the illumination tilt angle 



k x A. 



by multiplication of the incident wave function 
by a phase grating & m \x) = exp[iaV^ m \x)]) 
(cr = 2vr (m e c 2 + E ) / [XE (2m e c 2 + E )] being the 
electron-potential interaction constant) at the position 
of each slice and Fresnel propagation between the slices. 
We will assume an incident plane wave Aexp(27rifc • r) = 
Aexp[2Tri(k x x + k z z)\ arriving at the sample at an an- 
gle 9 — sin -1 (k x /k z ). Since holographic experiments can 
measure relative phases but not the absolute phase of an 
electron wave functions we will define z = in the plane 
of the exit face wave function, fixing the z-dependence 
of the wave function at the entrance surface. Combining 
this z-dependent phase factor with the scaling factor A 
we obtain (a;) — A^"') exp[2mk x x] at the entrance 
surface of the specimen. The electron wave function at 
the exit surface can then be obtained by the real-space 
multislice formalism Q as 
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Here A = \/k z — hc/y^Eo (2m e + Eq) is the electron 
wavelength, Eq the accelerating voltage, h Planck's con- 
stant, c the speed of light, m e the mass of an electron, 
A x the Laplace operator and the gradient in the re- 
direction. 

Implementing the multislice algorithm on a computer 
requires a discrete grid in the lateral dimension in addi- 
tion to the discrete slices in the z-direction. Keeping the 
lateral sampling distance S x constant we may simplify the 



notation by defining a dimensionless imaginary parame- 
ter 7fc x = 2mS x k x and replacing x by j x 8 x , V(x) by * Ja ., 
V^ n \x) by v£ m) , & m) (x) by and exp[27rifc x a;] by 

exp[7fc x jz]. This makes also the translation of the fol- 
lowing equations into a computer program a bit more 
apparent. 

Expanding the Fresnel propagation operator by using 
the relations 
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in orders of approximation of the Ewald sphere curvature 
given by the parameter f3 — ieX/ (AttS 2 ^. The sampling 
distance S x is typically quite a bit smaller than the image 
resolution defined by aberrations of the electron optics, 
the stability of the microscope, and the size of the ob- 
jective aperture (e.g. 0.1...0.5A for HRTEM at 1 A 
resolution). 

By chosing the appropriate real-space sampling, slice 
thickness and electron accelerating voltage one can force 
to have a fairly small value. Plugging the expansion 
of the Fresnel propagator (j2]) into {1} we can expand the 
expression for the exit face wave function in orders of 
0. The first order expansion neglecting, for the moment, 
terms of order 2 and higher looks as follows: 
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FIG. 2: (color online) Diagram of first (A) and second order 
(B,C) approximation of the Ewald sphere curvature effects in 
the real-space multislice algorithm for k x = 0. The A x op- 
erator is depicted by a node merging 3 paths, and the A x 
operator by a node merging 5 paths. Multiplying the scatter- 
ing potential along the the read solid, blue dashed and green 
dotted paths in A will produce all the possible monomials 
that are linear in j3, i.e. which involve only a single op- 
erator. Monomials resulting from paths shown in B and C 
(only one possible path is shown in each) are proportional to 
P 2 and involve either 2 successive A x operations, or a single 
A x operation. 



Note here, that the first term, the zeroth order expansion 
is independent of /3 and is just the well-known phase ob- 
ject approximation which includes multiple scattering to 
N th order but neglects effects due to curvature and tilt 
of the Ewald sphere. 

Figure^ illustrates expression graphically for an 
object partitioned into 4 layers. Examples of the 2 possi- 
ble scattering paths that scale as 1 are shown in figure 
GQj and [2b- While the conventional multislice algorithm 
[13| iterating between real and reciprocal space naturally 
includes all orders j3 n , in real-space algorithms they need 
to be included explicitly. Since increasing n slows down 
the computation enormously, terms 0(/3™- 3 ) in the com- 
putation of a single Fresnel propagation are usually ne- 
glected. An algorithm approximating the Fresnel prop- 
agation up to 0((3 n=1 ) will only have contributions of 
0(P n - 2 ) to the final wave function of the type shown in 
figure [5b, i.e. multiple operations, but neglect the 
equally important contributions shown in figure [2ji. 

The aim of this letter being the determination of the 
3-dimensional scattering potential from the observed im- 
ages we must consider the contributions to the exit face 
wave function in the order of significance, i.e. all terms 
up to a given order 0(/3"). Taking a closer look at figure 
[2] makes clear that in order to reconstruct N slices one 
needs to expand the Fresnel propagation up to at least 
n = (N — l)/2 (rounding up for even N). The result- 
ing system of polynomial equations of degree N is sparse 
and may be solved using globally convergent algorithms 
for solving multivariate polynomial sets of equations (see 



e.g. 14 1 with an example application given in [l5(). 

Holographic methods such as off-axis electron holog- 
raphy [llj , but also focal series reconstruction (see, for 
a very recent example, [l2j and references to 9 alterna- 
tive algorithms therein) are able to reconstruct the com- 
plex exit face wave function for each incident beam tilt. 
If, as is common practice, in the off-axis holographic 
reconstruction the side band is properly centered, and 
in the focal series reconstruction the images are aligned 
by cross-correlation or similar methods, neither of these 
methods will reconstruct the global phase shift exp^fc^ j x ] 
that is associated with the tilted illumination (the rela- 
tive phase factors in the second and third term of ([3|) 
remain, though). We can therefore drop this term allto- 
gether. However, a reference point common to wave func- 
tions of different incident beam tilt must still be chosen 
in order to fix the complex coeffcients A^ kx \ A vacuum 
region or another area of well-known scattering proper- 
ties within the field of view would be ideal, but if no 
such reference point can be defined the A^ k ^ parame- 
ters can also be included in the non-linear reconstruction 
algorithm. We will assume A^ k "' — 1 in the example 
presented below. 

For demonstration purposes, and in order to reveal the 
principles of the proposed methodology independent of 
the performance of a given polynomial equation solver 
we approximate the system of polynomial equations ([3|) 
by expanding the phase grating in terms of the po- 
tential assuming the validity of the weak phase object 
approximation. For a structure that has been split into 
N = 3 layers we approximate 
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converting (J3]) into a linear system of equations, the so- 
lution of which should be straight forward. It turns out 
the the matrix defining the resulting system of equations 
is slightly rank-deficient, because dropping all the non- 
linear terms in introduces 2 additional degrees of 
freedom, the overall phase relationship between the 3 lay- 
ers. Defining the relative potential of 1 column of pixels 
(e.g. in the example shown in figure [3] I defined 3 addi- 
tonal linear equation that force the potential to be zero 
(vacuum) for the right most pixel in each of the 3 layers) 
will remove this degree of freedom and produce, in the 
absence of noise, a unique reconstruction. The system 
of linear equations has been solved using singular value 
decomposition. 

Figure [3] a shows the phase shift of a trial structure 
that has been sliced into 3 equidistant slices. Scattering 
factors and length scales have been scaled to those of the 
gold (110) lattice, in order to mimic actual experimen- 
tal conditions. Since no noise has been included in this 
test, it is not surprising that reconstruction and original 
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FIG. 3: Phase shift cj> {m) (x) = crV (m) (x) for a 3-layered 
"phantom" structure (a) and the phase shift reconstructed 
from it using the linear approximation to © (b) demonstrat- 
ing a 3-dimensional reconstruction at atomic resolution. The 
scattering factor and lattice constant correspond to that of 
gold. The "defect" sites have been created by simply scaling 
the scattering factor of gold by a number less than 1 (indi- 
vidual defects from left to right: second row: 0.5, 0.5, 0, 0.7, 
third row: 0.3). Apsorption has been included by multiplying 
the potential V i - Tn \x) by (1 + 0.1*). Simulation parameters: 
object potential grid: 240 x 3 pixels, S x = 0.128A, e = 2.04A, 
E = 60kV, cr = 0.0011 (VA)~\ P = 0.49*, 5 different tilt an- 
gles of the incident beam: 9 = -1° (-y^ = -0.29*), 9 = -0.3° 
(j kx = -0.086*), 9 = 0° ( 7fc;r = 0i), 9 = 0.7° (j k:c = 0.20*), 
= 1° {lk m = 0.29*). 



look identical. The fairly low accelerating voltage of only 
60kV causes a (multiple scattering strength), (3 (Ewald 
sphere curvature) and ^jk x (illumination tilt sensitivity) 
to be of large enough values to make a 3-dimensional 
atomic resolution reconstruction of nano structures fea- 
sible, even in the presence of noise and for small beam 
tilt angles. 

Some iterative tomographic reconstruction algorithms 
allow the specimen geometry to be used as a constraint, 
in principle being able to reconstruct an object consist- 
ing of only N distinct layers to be reconstructed from N 
projections only. The angle between these N projections 
must, however, be quite large. Although requiring the 
same minimum number of tilt angles, in contrast to (lin- 
ear) tomography, where the 3-dimcnsional information is 
introduced by projecting along different directions, mak- 
ing use of dynamical scattering as proposed here, implies 
a tilt-angle dependent (complex) weighting factor for a 
set of monomials in the polynomial system of equations 
([3]). These weighting factors depend strongly on the ac- 
celerating voltage and can therefore be tuned to the scat- 
tering strength of the material, available beam tilt range 
and desired resolution in the direction parallel to the elec- 
tron beam. Modern TEMs being able to achieve atomic 
resolution imaging at Eq < 60kV (e.g. [T3|) will therefore 
be able to image the 3-dimensional potential distribution 
within nanostructures at atomic resolution without hav- 
ing to tilt the specimen. 



It should be noted that exit wave reconstruction for 
several illumination tilt angles is only one of several meth- 
ods for acquiring the experimental data required for the 
3D reconstruction described above. Alternative methods 
include Ronchigram focal scries recorded for different il- 
lumination or specimen shifts and, of course, holographic 
experiments at different specimen tilts. 

Summarizing, a method for reconstructing the 3- 
dimensional electrostatic potential of a TEM sample has 
been proposed. The method is based on a reformulation 
of the real-space multislice formalism for computing the 
multiple scattering of a fast electron within a TEM sam- 
ple in terms of a polynomial set of equations, identifying 
and keeping the most significant terms. The solution of 
the resulting sparse multivariate polynomial system of 
equations can be found by applying polynomial equation 
global optimization algorithms available in the literature. 

In the special case that the weak phase object ap- 
proximation is valid the polynomial system of equation 
transforms into a linear system of equations which can be 
solved by standard methods. Applying this simplification 
a 2-dimensional complex test object consisting of 3 layers 
has been reconstructed successfuly from 5 1-dimensional 
electron wave functions simulated for illumination tilt an- 
gles ranging from —1° to +1°. The construction of an 
efficient global optimization algorithm specialized for the 
particular type of systems of polynomial equations de- 
scribed in this letter is planned for the near future. 
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